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Abstract. In this article we propose and investigate a hierar¬ 
chy of mathematical models based on partial differential equations 
(PDE) and ordinary differential equations (ODE) for the simula¬ 
tion of the biophysical phenomena occurring in the electrolyte fluid 
that connects a biological component (a single cell or a system of 
cells) and a solid-state device (a single silicon transistor or an ar¬ 
ray of transistors). The three members of the hierarchy, ordered by 
decreasing complexity, are: (?) a 3D Poisson-Nernst-Planck (PNP) 
PDE system for ion concentrations and electric potential; {ii) a 2D 
reduced PNP system for the same dependent variables as in (i); 
(in) a 2D area-contact PDE system for electric potential coupled 
with a system of ODEs for ion concentrations. The backward Eu¬ 
ler method is adopted for temporal semi-discretization and a fixed- 
point iteration based on Gummel’s map is used to decouple system 
equations. Spatial discretization is performed using piecewise lin¬ 
ear triangular finite elements stabilized via edge-based exponential 
fitting. Extensively conducted simulation results are in excellent 
agreement with existing analytical solutions of the PNP problem 
in radial coordinates and experimental and simulated data using 
simplified lumped parameter models. 


Keywords: Bio-hybrid systems; neuro-electronic interfaces; multi¬ 
scale models; electrodiffusion of ions; functional iterations; numerical 
simulation; exponentially fitted finite elements. 


1. Introduction and Motivation 


In this article we address the study of a class of problems arising 
in the context of Bioelectronics, a recently emerged discipline at the 
crossroad among Nanotechnology, Solid-State Electronics, Biology and 
Neuroscience. The focus of our investigation is on the mathematical 
and computational modeling of bioelectronic interfaces (see IHKITj for 
a review and [26l [20], HG] |3l] for a selection of signihcant applications). 
Bioelectronic interfaces are bio-hybrid structures constituted by living 
cells attached to an electronic substrate and surrounded by an elec¬ 
trolyte bath. An example can be seen in Fig. 1(a) which shows an 


electronmicrograph of hippocampal neuron cultured on an electrolyte- 
oxide-silicon held-effect transistor (EOSFET) |12]. Fig. 1(b) reports 
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a schematic cross-section view of a neuro-chip, which allows to iden¬ 
tify the main parts of the bio-hybrid system: the cell, the extracellular 
bath, the thin interstitial cleft separating the cell and the electronic 
substrate, the protective oxide layer deposited on the top of the sub¬ 
strate, and the source-to-drain transistor structure. 



(a) Rat neuron on electronic substrate 



(b) Neuro-chip 


Figure 1 - Left: rat neuron grown on an EOSFET, image 
reprinted from |42j . Right: schematics of a neuro-chip, image 
reprinted from |35j . 

In the basic function mode of the EOSFET, as a consequence of 
the cellular activity elicited by the application of an external stimulus, 
ionic current flows through the adhering cell membrane and along the 
cleft. The resulting extracellular voltage turns out to play the role of 
the gate voltage which controls electrical charge flow in the substrate 
and, ultimately, the current flowing out from the drain terminal of the 
device. 

The interface contact in the scheme illustrated in Fig. [T] is realized by 
the thin conductive electrolyte separating the two subsystems, whose 
amplitude is smaller than the cell radius by about three orders of mag¬ 
nitude. Therefore, the cell-chip junction forms a planar electrical core¬ 
coat conductor and the main physical phenomena (ion electrodiffusion 
and gate voltage modulation) take place in this three-dimensional re¬ 
gion whose vertical thickness is much smaller than the two-dimensional 
area where the cell adheres to the substrate. 

The above description indicates that a sound mathematical picture 
of a bioelectronic interface requires the adoption of a genuine mul¬ 
tiscale perspective. For this reason, in this article we continue our 
analysis started off in [Ij and numerically investigated in [5] and pro¬ 
pose a hierarchy of models based on PDFs and ODEs for the sim¬ 
ulation of the biophysical phenomena occurring in the 3D interface 
contact described above. The hierarchy includes the following three 
members, ordered by decreasing level of complexity: {i) a 3D Poisson- 
Nernst-Planck (PNP) PDF system for ion electrodiffusion and electric 
potential dynamics |39]; (m) a 2D reduced PNP system for the same 
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dependent variables and phenomena as in (z); {in) a 2D area-contact 
PDE system for electric potential dynamics coupled with a system of 
ODEs for ion dynamics. This last member of the hierarchy is a variant 
of the area-contact model proposed and studied in [S]. 

Model (i) is the most accurate in the hierarchy but, of course, re¬ 
quires a considerable amount of computational effort for its numerical 
simulation. Model (m) is obtained by averaging the 3D PNP equations 
in the direction z perpendicular to the electrolyte cleft. This model re¬ 
duction procedure leads to a modihed PNP system to be solved in a 2D 
plane x-y parallel to the substrate. Model {in) is a further reduction of 
{a) obtained by neglecting spatial dependence of ion concentrations in 
the electrolyte cleft. This leads to a time-dependent 2D Poisson equa¬ 
tion for electric potential coupled with electrolyte cleft ion dynamics 
described by a system of ODEs as done in [S]. In all members of the 
hierarchy, iono-electric coupling between substrate and electrolyte is 
accounted for by “lumped” transmission conditions expressing conti¬ 
nuity of dielectric and ionic fluxes across the interfaces. Electrodiffusive 
ionic coupling between cell(s) and electrolyte is described through a va¬ 
riety of transmembrane currents including the Goldman-Hodgkin-Katz 
and Hodgkin-Huxley models [32l [M] . 

The backward Euler method is adopted for temporal semi-discretization 
and a hxed-point iteration based on Gummel’s map [23 is used to de¬ 
couple system equations. Spatial discretization is performed using a 
generalization to axisymmetric cylindrical coordinates of the piecewise 
linear triangular hnite element scheme stabilized via edge-based expo¬ 
nential htting proposed and analyzed in |46j . 

Extensively conducted simulations using the full 3D PNP model re¬ 
veal that ion concentration and electric held variations mainly occur 
in the electrolyte cleft. Sensible results are also obtained addressing 
non ideal effects occurring in realistic devices, such as undesired cel¬ 
lular activity detection on more than one electrode and cellular cross¬ 
stimulation. Simulation experiments using the 2D formulations pro¬ 
posed in the present article demonstrate their excellent agreement with 
experimental and numerical results in the existing literature [361 [8] but 
with a signihcant saving of computational effort with respect to the 
solution of the full 3D PNP system. 

A short outline of the article is as follows. In Sects. |2] and [3] we 
illustrate the hierarchy of mathematical models used to represent ion 
electrodiffusion and electric potential distribution in the bioelectronic 
structure. In Sect. we describe the numerical techniques adopted 
for the solution of the discrete problem. In Sect. we address the 
validation of the proposed computational model in the simulation of 
several test cases of biophysical signihcance. In Sect, [^we summarize 
the main contents of our analysis and indicate some perspectives for 
future research directions. 
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2. Three-Dimensional Model 

In this section we illustrate a three-dimensional model of ion elec¬ 
trodiffusion throughout the interstitial cleft separating the cell and the 
electronic substrate under the application of an external stimulus. 


2.1. Geometrical model. Fig. 2(a) shows a 3D schematic picture of a 
cell-to-substrate interface. The cell shape is represented as a rotational 
solid separated by the planar substrate by a thin electrolyte domain 
Despite the extracellular fluid is surrounding all the cell, most of the 
effects resulting from cell stimulation occur in the adhesion region as 


is demonstrated in the numerical experiments reported in Section 5.2 


For this reason, we restrict the geometrical space where to solve the 
mathematical problem to the 3D region Dez, see Fig. 2(b), which in¬ 
cludes the cleft between the attached cell membrane and the device, but 
also the part of electrolyte in the neighborhood of the cell close to the 
substrate. The layer thickness 5j is of the order of 50 -P 100 nm while 
the cell radius is around 10 pm in the considered applications [SI IHU]. 




(a) Cell, electrolyte and substrate (b) 3D computational domain and inter¬ 
faces 

Figure 2 - Left: geometrical model of the whole bio-hybrid sys¬ 
tem. A cell surrounded by an electrolyte bath is attached to an 
electronic device. Right: the computational domain is the 
thin layer of electrolyte between the cell and the substrate. Seven 
different boundary regions are distinguished: the upper surface is 
divided in Tceii (the cell attachment area) and T^f (surface sepa¬ 
rating the electrolyte cleft from the extracellular fluid). 
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2.2. The Poisson-Nernst-Planck system. The Poisson-Nernst-Planck 


system (PNP) for ion electrodihusion reads [39]: 


(la) 

uC' 

-h divfj {ci,ip) = 0 

i = 1,..., M 

(lb) 

fi (ci, ip) = -DiVci /ij-A-CjE 

\^i \ 

i = 1,..., M 

(Ic) 

D, = 

Gil 

i = 1,..., M 

(Id) 

divE = -p 
e 


(le) 

E = —V(p 

M 


(If) 

p = q'^ZiCi. 
i=l 



Eq. (la) is the continuity equation describing mass conservation for 
each ion whose concentration is denoted by Cj i = 

M > 1 being the number of ions flowing in the electrolyte fluid. Each 
ion flux density f* (m“^ s“^) is dehned by the Nernst-Planck rela¬ 
tion (lb) in which it is possible to recognize a chemical contribution 
and an electric contribution, in such a way that the model be regarded 
as an extension of Pick’s law of diffusion to the case where the diffusing 
particles are also moved by electrostatic forces with respect to the fluid. 
The quantity Zi is the valence of the Tth ion, while fii and Di are the 
mobility and diffusivity of the chemical species, respectively, related by 
the Einstein relation (Ic) where Vth = ksT/q (V) is the thermal poten¬ 
tial (fcs is the Boltzmann constant, T is the absolute temperature and 
q is the elementary charge). The electric held E (V m“^) due to space 
charge distribution p in the electrolyte is determined by the Poisson 
equation (Id) which represents Gauss’ law in differential form, e being 
uthe dielectric permittivity of the huid medium. For further analysis, 
it is useful to introduce the electrical current density j* (A m“^), equal 
to the number of ion charges howing through a given surface area per 
unit time and dehned as 


(Ig) ji := qzi^i = -qZiDiVci qpi\zi\ciE i = 1,..., M. 

Remark 1. The PNP system ([^ has the same format and structure as 
the Drift-Diffusion equations for semiconductors (see, e.g., ^\), but 
it is applied to a different medium (water instead of a semiconductor 
crystal lattice) and includes, in general, more charge carriers than just 
holes and electrons, as in the case of semiconductor device theory. 


2.3. Boundary and initial conditions. Let t and x denote the time 
variable and the spatial coordinate, respectively. We denote also by 
P := dVLei the boundary of the computational domain in Fig. 2(b) and 
by n the unit outward normal vector on P. 
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The initial conditions c°(x) = Ci(0,x) and V5°(x) = (^(0,x) are de¬ 
termined by solving the static version of the PNP system ([^ in the 


domain hlgo which corresponds to setting ^ = 0 in (|la|) for each ion 


i = 1,..., M. 

The boundary conditions deserve a deeper discussion because they 
need to mathematically express the coupling between the electrolyte 
cleft and the surrounding environment, comprising the extracellular 
fluid and the active parts of the bio-hybrid system, namely, the cell. 


the membrane and the electronic substrate (see Fig. 2(a)). Referring to 
Fig. 2(b) I for the notation, we distinguish among seven different regions 


in the boundary F: four lateral sides Fj, i = 1,...,4, a lower face F^^b 
in contact with the substrate, and an upper face divided into two parts, 
the area attached to the cell Tceiu and the free region Fe/, covered on 


top by the surrounding volume of extracellular huid. Accordingly, the 
following conditions are enforced on the electric held and the particle 
huxes: 

(2a) 

f — Vbath 

onFiUF2UF3UF4 

(2b) 

[D ■ nlr., = 0 

onTef 

(2c) 

[D ■ n|r„„ = 0 

on F cell 

(2d) 

[D . nlr.., = 0 

on 

(2e) 

^ _ ^bath 

onFiUF2UF3UF4 

(2t) 

|fi ■ n]r,^ = 0 

onTef 

(2g) 

[fi ■ = 0 

on F cell 

(2h) 

fi ■ n = 0 

ou 


having denoted by 1-]^ the jump operator restricted to the interface C. 


Eqns. (2a) and (2e) are Dirichlet boundary conditions that can be 


interpreted as “far held conditions”, meaning that sufficiently far from 
the surface where the cell is attached to the substrate, we can assume 
the electric potential (p to be hxed at a given constant reference value 
and each ion concentration to be hxed at a given constant value 
The quantities i = 1,..., M, are physiologically given values in 
such a way that the electrolyte is a neutral solution, i.e.. 


M 


(2i) 


P Pbath Q ^ ^ 


Mth ^ Q_ 


2=1 


Cleft-cell coupling. The most relevant physiological phenomena occur¬ 
ring in the bio-hybrid interface depend on the properties of the cell 
membrane and can be modeled as the sum of two contributions, one 
from the lipid portion of the membrane and one from the ionic channels. 
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(b) Lumping in Tceil 


Figure 3 - On the left: cell and electrolyte separated by the 
membrane with its physical thickness tM- On the right: cell and 
electrolyte separated by an interface Tceii with zero thickness, re¬ 
sult of the lumping of the original boundaries Fi and r 2 of the 
membrane region. 


The membrane subdomain (shown in Fig. 3(a)) has a thickness Im 
of the order of 5 -F 10 nm, which is much smaller than the characteris¬ 
tic size of the domain VLei- Therefore, a geometrical discretization of 
this small region may give rise to a huge number of degrees of free¬ 
dom for the numerical method. To reduce computational complexity, 
we apply the membrane model proposed in [32l E]. This amounts to 
assuming that ip varies linearly across the membrane thickness so that, 
upon introducing the two dimensional manifold Tceii corresponding to 
the middle cross-section of the membrane volume, the transmission 


condition (2c) across the two dimensional manifold Tceii becomes 


(3a) Dc ■ ric = -De ■ Ue = -em 


^m2 ^ml 


t 


M 


— Cm i^Pm2 Ccell) , 


where Em is the membrane permittivity, Vceii is the intracellular po¬ 
tential while (prni and ipm 2 are the traces of p at both sides of Tceii 


(cell and electrolyte, respectively). Condition (3a) expresses a capac¬ 
itive coupling between cell and cleft through the membrane specific 
capacitance Cm '■= ^m/^m (Fm“^). 

To account for membrane channels, the ionic current in the case of 
active cells is described by the following generalized Hodgkin-Huxley 
(HH) model (for a detailed description see 


(3b) 




where s is a vector collecting the gating variables n, m and h responsi¬ 
ble of channel opening probabilities, while and c are arrays of size 
M containing all the ion concentrations inside and outside the cell. If 
instead only passive cells are considered, as in most simulations per¬ 
formed in this work, the adopted model is the Goldman-Hodgkin-Katz 
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(GHK) equation for the current density 
(3c) 


Ji”" = PiZiq 


Be 


iy^cell P) 


„cell 


-Be 


Zi {Vcell - (f) 


Vth J ^ \ Vth 

where Pi is the permeability constant of the specihc ion (ms“^) and 
Be{x) := x/{e^ — 1) is the inverse of the Bernoulli function. With 
these descriptions of the transmembrane current density j*”* for each 
ion species, condition becomes 


(3d) 


tm 


n, 


= -f* ■ n. = - — 


k 

qz. 


Cleft-substrate coupling. In the present work, the action of the electro¬ 
chemical bounding of ions at the interface between electrolyte cleft 
and substrate is neglected, and the semiconductor device is assumed 
to behave as a MOS capacitor mathematically described through a 
lumped equivalent model as in (3a). Referring to Fig. |^the capacitive 
coupling on Fs^f, is 


(4a) 


D, ■ n, = -Dp ■ rip ~ -e 


— Psl 


— —Cs {^32 — Vq) , 


where and ips 2 sue the traces of ip on both sides of The 

function Vq = Vg (t) denotes the value of the potential on the gate 
contact, taken to be spatially constant according to the hypothesis 
of ideal metallic behavior of the gate. Regarding the particle fluxes. 


condition (2h) already states that there is no current injection in the 


electronic device. 



Figure 4 - Electronic substrate model and coupling with the electrolyte. 


Electrolyte-electrolyte artificial coupling. As shown in Fig. 2(a), the 


electrolyte domain is restricted to a thin sheet of amplitude 5j. This 
approximation leads to the boundary conditions (2b) and ([^ on the 
hctitious boundary Fp/. Assuming again that the potential is a linear 
function of z, we can rewrite (2b) as 

(4b) Dgxt ■ ^ext Dint ' ^int — C {}Pint Pext) — C {jpint ^bath) i 
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where (fint and tpext are the traces of tf respectively on the two sides 
of he/ (see Fig. [^. The reduced model (4b) consists of assuming that 
far away from the boundary Fg/ the potential is at the reference value 
hfeai/i, C* being a hctitious capacitance introduced to relate the value of 
the potential in the electrolyte, inside and outside of the computational 
domain A possible modeling approach to estimate the value of C* 
consists in taking a fraction 1/k oi the value of Cm- Computational 
experiments indicate that k = 5 is an appropriate choice. 



illustrating the coupling condition between tleZ and the external 
remaining electrolyte enforced on T^f. 


Regarding particle fluxes, we assume that, far away from Fgj, ion 
concentrations can be considered to be equal to their bath value 


and the electrolyte to be electroneutral. Then condition (2f) can be 
rewritten as 
(4c) 


^ext ^ext 


^int'^int 


= -v* (cf* - ^ V* (cr* - ^ = 1, 


,M. 


The above relation is a Robin condition for the particle flux density 
which physically expresses the fact that ions are allowed to cross the 
hctitious interface Fg/, as it should be in the non-truncated electrolyte 
domain. Mathematically, the quantity v* is an effective permeabil¬ 


ity (ms whose value can be estimated by equating hux (4c) to a 


fraction 1/k* of the hux through the membrane (3d). Computational 


experiments indicate that k* = 20 is an appropriate choice. 


3. Hierarchical models 

The analysis of [8l |36] shows that the ion current density entering 
the cleft through the membrane mainly hows parallel to the 2 ;-axis. 
Once inside the cleft, the direction of the current density changes into 
the radial one. The time needed by the ions to how across the cleft 
thickness is of the order of 10“^ s. Particles move many times up and 
down along the ^-direction because the ratio between the cleft thickness 
and the radius of the attached area is of the order of 10“^ and the 
resulting contribution of this random motion to the vertical current is 
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equal to zero. Thus, it appears to be reasonable to derive a family of 
two-dimensional models in the x — y plane from the 3D PNP system 
of Sect. 1^ This is the object of the discussion below. 

3.1. Model reduction: from fully 3D to 2.5D ion electrodif¬ 
fusion . We place a coordinate system with the origin in the middle 
of flei- The plane in the middle of the cell-chip junction, depicted 
in Fig. is going to be the new two dimensional domain Q 2 D'- it is 
equidistant from Ps^f, and from TefUTceii, which are respectively placed 
at z = —Sj/2 and z = +Sj/2. 



Figure 6 - Schematics for the geometrical reduction in the x-y 
plane: the middle plan of the cleft becomes the two-dimensional 
domain 0,20 ■ Vxy is a control volume used to compute the integrals 
and the fluxes. 

We follow a procedure based on the integration of the three-dimensional 
Eqs. Q on the test volume of Fig. This latter is a parallelepiped 
of volume Vxy = Sjhxhy, where hx and hy are inhnitesimally small. 
Introducing the following integral means: 


(5a) 



and integrating the continuity Eq. (la), we obtain 


(5b: 
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fa;, iy and being the components of the fluxes in the three directions. 
An analogous result is obtained for the Poisson equation (Id). 

The values of the fluxes f* and of the electric displacement D on 
and Sfooi are unknown and need be computed according to the boundary- 
conditions applied on these surfaces. The boundary conditions on 


and ^bot 

can 

be written as: 





(6c) 


■ n = 

fr {t\ 

top 


■ 9^) 

on 

^top 

(5d) 


■ n = 

t‘ (*; 

hot 


^bot') 9^) 

on 

^bot 

(5e) 

D 

■ n = 

gtop (fi 



9^top; 9^) 

on 

^top 

(5f) 

D 

■ n = 

9bot (L 

bot 


9^foot: 9^) 

on 

^bot: 


where the functions /f"*, gtop and gi,ot depend on the “averaged” 
quantities defined in (5a), but also on the quantities evaluated on the 
surfaces S^p and T,hot, defined as: 

(5g) 

(5h) 


top 
a- := Ci 


I Stop 

^top ■= 


:= c 


top 


* I ^bot 

^bot '■= 


Dividing (5b) and the analogue for the Poisson equation by |14p| = 
5jhxhy, and taking the limit as hy —)■ 0, we obtain the following 
averaged 2.5D PNP model in Vt 2 D'- 


(6a) 

(6b) 

(6c) 

(6d) 


— + div f • + — = 0 


f,; = -A V 


xyCi T ~ 

Vth 


— 1 1 ^ _ 

diVa;pD T ~gtop T ~^9bot Q / ^ 

3 J 2 

D = -eVxy^- 

The boundary conditions applied on dfl 2 D simply reduce to: 

(6e) Q = 

(6f) 9^ ^bath’ 

The reason for qualifying the novel electrodiffusive model @ as a 
2(+l/2)D=2.5D formulation is related to the definition of the source 
flux terms and gtop^dbot: object of the next section. 


3.1.1. The reason for the ”2.5D”: boundary layer model approximation. 
Let xc, denote the characteristic function of a domain C Then, 
since the upper surface of Vt^i is the union of two different parts, Vceii 
and Pe/ (cf. Fig. [^, we decompose and gtop into the sum of two 
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Figure 7 - Cross section in the x-z plane of the three- 
dimensional cleft, showing the scheme for the model of 
the boundary layers near the surfaces T^eii and Fs^f, 
{H is the layer amplitude). Decomposition into three 
subdomains: = {{x, z) s.i. z ^ [5j/2 — H] 5j/2]}, 

^2 = {(x, z)s.t. z &[—5j/2 + H] 5j/2 — H]} and 

D 3 = {(x, z) s.t. z G [—5j/2-, —5jl2 + H]} 


contributions as: 

(7a) C = + nio^XT., 

(7b) Slop = +9)4,,\r.,. 

The coupling conditions on Vceii and T^ub are between different en¬ 
vironments and typically give rise to the occurrence of boundary lay¬ 
ers |H [TSl 133]. These latter are in the form of electrical double layers, 
of which we only account for the diffuse layer, neglecting the ions at¬ 
tached to the surfaces, as in the Gouy-Chapman approximation cai. 
In Fig. [3 we focus our attention on a x-z cross section of the whole 
three-dimensional electrolyte cleft aX y = y, denoted VL^z- This latter 
is partitioned into three distinct subdomains as Vt^z = U f22 U fls, 
H being the amplitude of the two boundary layer regions f22 ami 
According to physical evidence, for every fixed point x of the x axis, 
we assume that 


dip{x, z) 
dz 


dci{x,z) 

dz 


0 infl2. 


and we set ip{x,z) = Tp{x,z) axidciix^z) = Ci{x,z) for all 2 ; G ^ 2 - 
These two definitions amount to extending along the zi-direction (in 
the sole interval ^ 2 ) the averaged values determined by the model de¬ 
scribed in Sect. 113 We also introduce further assumptions on the 
electric potential and the particle fluxes in the two boundary layer 
subdomains: 

(1) (p is linear in fli and fls and continuous at z = 5j/2 — H and 
at z = —Sjl2 -b H] 

(2) fj is constant in fli and fls. 
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Figure 8 - Scheme of the distributions in the ^-direction at a fixed 
X for z) and for Ci{x, z) (where the considered ion is positively 
charged) for the cross section depicted in Fig. [tI 


The spatial distribution of <p(x, z) for a fixed point x is schematically 
depicted in Fig. 8 (a) Assumption 1. indicates that the electric held 
is piecewise constant over VL^z (and equal to zero in ff 2 )- Also the 
particle huxes are piecewise constant over Vt^z (and equal to zero in 
because both drift and diffusion terms are null there). In order to 
determine the concentration Ci{x,z), we integrate the Nernst-Planck 
transport equation (lb) in ffi and fls- The resulting distribution of 
ions is piecewise exponential over VL^z, continuous at z = 5j/2 — H and 
at z = —5jl2 + H, and constant in fl 2 , as depicted in Fig. 8 (b)[ The 


corresponding mathematical expressions for the boundary huxes 
and fbot are: 




Zii^Ptop top Tj f top _ 

Wh ; “ V 


<P- 


c“ - Be (- 

Vth J \ 


Zi{ip - ipbot) \ _ 
Cj 


Vb 


th 


The above described modeling reduction procedure is equivalent to 
applying the Scharfetter-Gummel (SG) exponentially htted approxi¬ 
mation in Gi and G 3 001 . Regarding the electric displacement, with 
the above approximation we obtain: 


(9a) 

(9b) 


cell _ ftx)p__^ 

ytop ^ TT 


Qbot c 


H 


To determine the values if top and ipbot we use ([^ into the 

coupling boundary conditions (3d) and ( 2 h), and (|^ into the coupling 
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boundary conditions (3a) and (4a), to obtain: 


(10a) 


top\ 

S r 77 

1 

Be ( 

(10b) 

„bot 1 _ 

* 

(10c) 


(lOd) 



Be(- 


CiBe {zi{(ptop - ^)/Vth) + 


Zi{}P ^bot)/yth) — 
Ci 


jT^\ 

qziDi) 


c 


Be {zi{ip - ipbot)/Vth) 
CmVccU + 


M -\- ^/H 
1 

C's + e/H 


{CsVg + . 


Let us now consider the artificial surface Tej. No boundary layer is 
expected to occur there, so that we simply set 


( 11 ) 




I _ — 


- e/ 


= Ci 


and use these values into the electrolyte-electrolyte coupling condi¬ 
tions ( [4^ and (4b) to compute the functions and f^J^p introduced 
in ([^. The final 2.5D ion electrodiffusion model in ^20 then reads: 

L"“) Xir., = 0 


(12a) 

f)r -1 11 

+ div./, + -/‘2„ xir„„ + j:ft“ + T'’" 

(12b) 

(12c) 


f-i xy^i -rr xy^') 

yth 


div^^D + ^gZp xli 
Cj 


(12d) 


+ YStot+ fc’ {<P - Uoii.) xir,, = 

Cj . 

D = -eVxyW- 


System (12) is completed by the same kind of initial conditions as for 
the 3D PNP system (^, by the boundary conditions (6e)-(6f) and by 
the set of relations (|8)-(11) 

3.2. A 2D electrical model for ion transport: the Area Contact 
formulation. In the same spirit as done for single cells by Fromherz et 
al. in [SI 133 ] and for multi-electrode arrays (MEAs) in [28], it is possible 
to obtain from the 2.5D equation system ([^ a genuine 2D ion transport 
model, denoted Area Contact model. To this purpose, following the 
approach of [S], we consider only the attached area as computational 
domain Vt 2 D-i we neglect the variation of (p and c* in the 2 ;-direction and 
we also assume that ion densities are spatially homogeneous. In this 
manner, the functions f^°^, gtop and gbot are still computed using 
the 3D boundary conditions of Sect. 2.3 but setting ptop = Pbot = P 
and = q. 
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Omitting from now the notation (•), we sum the M continuity equa¬ 


tions (6a) and use the Poisson equation and the boundary coupling 


conditions dehned in (3a) and (4a) to express the time derivative of p 
as 

- {p6j) = - (div,,D) + - {Cm^ + Cs^) - - {CMVcell + CsVg) . 

Replacing the previous relation into the sum of the continuity equations 
we end up with the 2D Area Contact model for ion transport: 

(13a) 


{Cm + Cs)^ + div,„ ( 




d 


dt J " jS + {CmV„u + CsVc) 


(13b) 


•cond _ 

otot 


\zi\ piCidjVp. 


For given ion concentrations q, system ( |13[ ) is a parabolic boundary 
value problem for the dependent variable p = p{x, y, t). Following the 
terminology adopted in |8] , we refer to system ( |I^ as the “2D electrical 
model”. Ion dynamics should also be accounted for, therefore at each 


time level we hrst solve (13) and compute the integral mean of p over 

^2D 


(13c) 


VAt) ■= 


fn,o f V’ 


|fi 


2D\ 


|D 2 _d| denoting the area of the contact area. Then, we use Vj{t) as an 
input voltage to solve the system of ODEs corresponding to the electri¬ 
cal equivalent circuit proposed in |S], to determine the concentrations 
Ci that must be passed to (13b) to advance to the next time level. 


Following again the terminology adopted in |H], we refer to this latter 
PDE-ODE coupled system as the “2D electrodiffusion model”. The 
physical accuracy of the 2D models introduced above are investigated 
in Sect. 15.6.31 


4. Functional Iterations and Numerical Discretization 

In this section we describe the techniques used to numerically solve 
the mathematical models introduced in Sections and The adopted 
strategy is composed of three steps: 

(1) Temporal discretization 

(2) Linearization 

(3) Spatial discretization 

Step 1 . For the temporal discretization we adopt the Backward-Euler 
scheme to approximate all time derivatives. Since in most applications 
considered in this work the input signal (usually the intracellular po¬ 
tential Vceiiit)) is expressed as a combination of Heaviside functions, 
the time step of temporal advancement At is a-priori appropriately 
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chosen in numerical simulations according to the following strategy: 
in correspondence of the on/off switching time of the signal, At is set 
equal to a small value, in the order of 1 x 10“® s; then, once transients 
are exhausted. At is suitably increased up to a value in the order of 
1 X 10-3 s. 

Step 2 . In order to handle the intrinsic nonlinearity of the models, 
we apply a functional iteration procedure that is widely used in the 
decoupled solution of the DD semiconductor device equations. The 
method is the well known Gummel Map [2Z1 EHl 1221 E], a staggered 
algorithm where each variable of the problem and its corresponding 
equation are treated in sequence until convergence. Gummel’s map 
in the steady state regime is theoretically investigated in the seminal 
book IZ7I to which we refer for details and further bibliography. Un¬ 
der proper conditions on problem geometry and boundary data, the 
Gummel decoupled iteration is proved to admit a hxed point and also 
to be a contraction. This implies uniqueness of the solution of the 
nonlinear PDE equation system under investigation. The convergence 
rate predicted by the analysis of [22] is linear, but computational ex¬ 
perience reveals that the Gummel algorithm in stationary conditions is 
exceptionally rapid and robust with respect to the choice of the initial 
guess, which makes it much preferable than Newton’s method despite 
its theoretically predicted lower order of convergence. 

Step 3. For the spatial discretization of the linearized PDFs we adopt 
the piecewise linear conforming Galerkin-hnite element method (G- 
FFM) stabilized by means of an exponential htting technique (Fdge 
Averaged Finite Element method (FAFF)) [21 HOl UHl H6] . in order 
to deal with possibly dominating drift terms and avoid the onset of 
spurious oscillations in the computed solutions. 

4.1. The EAFE method in axisymmetric geometries. In this 
section we illustrate the extension of the EAFE method proposed in [46] 
to treat the case of three-dimensional electrodiffusive problems in ax¬ 
isymmetric geometries as in the numerical simulations reported in Sect.]^ 


4.1.1. The electrodiffusion model problem. Let Gas = (0, i?) x (0,Z) be 
the computational domain shown in Fig. j^with boundary P := dQas 
divided into two disjoint subsets P d and P n and with outward unit 
normal vector n = [ur, Let u be the dependent variable (having 

the physical meaning of ion density), related to the ancillary dependent 
variable n through the constitutive equation 

(14a) u = ne^ 


where is a given function representing a normalized electric potential. 
The change of variable (14a) allows to write the DD flux density (lb) 
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symmetryl 


axis 



J 


Figure 9 - Schematics of an axisymmetric configuration. 


in an equivalent diffusive form and n undergoes the name of Slotboom 
variable [H]. The model of electrodiffusion of n that we consider in the 
present section is the following boundary value problem in self-adjoint 
form: 



(14b) 


where J(n) = [Jr{n) Jz{n)Y = iae^[dn/dr, dn/dzY, /i G C°(r2as) being 
the ion mobility such that /i = /i (t) > /iq > 0 Vx G hlas, and where -0 
is a continuous piecewise linear function over such that the drift 
held is b := V0 (a normalized electric held). We also assume the 
reaction coefficient c G L°°{flas), with c > 0 a.e. in the source 
term f E (flas) and the boundary hux jjy ^ (Tat). For any given 
functions 0 and 9 belonging to let us endow L‘^{flas) with the 

weighted scalar product (•, ■)^ 



(14c) 


where the weight function is cu = ^/r, cf) := cjcj), 9 : = u9 and {•, •) is the 
usual scalar product in L‘^{Qas)- Finally, let 


(14d) 4" := HY i^as) = {veH^ (nas) : = O} , 

endowed with the equivalent norm, using Poincare-Friedrichs inequality 
(see [38], Chapt. 1) 

(14e) ||«.||^.:=||Vm|L = ({Vm,V«.)J'''2 


Vw G V. 
















IMANUELA ABBATES MATTEO PORRO^’'^, THIERRY NIEUS^ AND RICCARDO SACCO^ 


Then, the weak formulation of problem (14b) reads: 
hnd n & V such that: 


(14f) 

where: 

(14g) 

(14h) 


(n, v) = (n) 


'iveV 


(n, v) = — (/re^Vn, Vn)^ + (^ce^n, 
Fuj (n) = (/, v)^ + / jNvds^ 

Jrff 


and where = rds is the curvilinear abscissa in radial coordinates, 
ds being the usual curvilinear abscissa. Using the Lax-Milgram lemma 
(cf. [3H], Chapt. 5) we can prove existence and uniqueness of the solu¬ 
tion n G U of ( 14f| ) and the following stability estimate [1] 

ll/llL2(n,,) + Cn IIjaI' 


\n 


\v< 


lL2(rjv) 






where Cp and are the Poincare’s and trace constants, respectively, 
while ipm is the minimum of '0 in f^as- Existence and uniqueness of n 
obviously imply the existence and uniqueness of the function u that 


solves the boundary value problem corresponding to (14b) but written 


as the continuity equation (la). 


4.1.2. Extension of the EAFE method to axisymmetric geometries. Let 
7h be a regular triangulation of the domain made by triangles K 
(cf. [3H], Chapt. 3). On Th we introduce the hnite dimensional subspace 
14 C U made of piecewise linear conforming hnite elements vanishing 
on r^). Denoting by G 14 the hnite element approximation of n, the 


application of the EAFE method to problem (14f) consists of hnding 
nh G 14 such that: 


(14i) {Ph, Vh) = {vh) 


\/Vh G 14 


where is an approximate bilinear form constructed in such a 

way that the approximation 3h{nh) of the hux J(n) over each trian¬ 
gle iP is a constant vector whose tangential component over each edge 
e G dK is computed by replacing the dihusion coefficient with its 

harmonic average along e. After computing the local stihness matrix 
as in the case of Cartesian orthogonal coordinates, A^ in still mul¬ 
tiplied by the weighting factor rdrdz (see [iiiiniE]). The following 
formula is used for the approximate evaluation of the local stihness 
matrix of the EAFE method in axisymmetric geometries 


(14j) 


A 


K 


fi 

0 

0 


0 

^2 

0 


O' 

0 

^3. 




^KeTh 


Ti being the radial distance of the midpoint of edge e* from the origin 
of the coordinate system, i = 1,2,3. The above approach proves to 
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be numerically stable and accurate as demonstrated in all the com¬ 
putational experiments reported in Sect. For the discretization of 
the local reaction and source terms, we adopt the same trapezoidal 
quadrature used in the Cartesian case and obtain 

0 0 


(14k) 


= 


'cie'^^ri 

0 

0 


C2e'^^r2 


0 

C3e^^r3_ 


3 


(141) F5 = 

/Ri 

Al + I 

O ^ O 

(iwg 

+ iiV,3)Fl|ei|4i 

f2r2 

(JY,3 

+ jAr,l)F2 62 5e2 



o Z 

{JN,! 

+ jAr,2)F3 63 (5e3_ 


VK e Th 


where r* is the radial distance of node i from the origin of the coordinate 
system, i = 1,2,3, \K\ is the area of K, \ei\ is the length of edge e*, 
i = 1,2, 3, while 5^^ is equal to 1 if G Fat and 0 otherwise, and, 
hnally, q, fi, jN,i and are the values of the Pi-interpolants of c, /, -0 
and jisf at each node i of K, i = 1,2,3 


Upon assembling (14j), (14k) 
and ( |141[ ) over the grid and applying the inverse of (14a) at each node 
of Th [SIE], we end up with the following linear algebraic system 

(14m) = F^, 

where u is the vector of nodal values of the dependent variable u while 
and Ftj are the global stiffness matrix and load vector of the EAFE 
method, respectively. As proved in pp, is an irreducible M-matrix 
with respect to its columns. This implies that (14m) admits a unique 


solution and that the discrete maximum principle holds under the same 
conditions as in the Cartesian case studied in [16]. In particular, if 
Fij > 0 (in the componentwise sense), then the solution of (14m) is 
such that u > 0. 


5. Numerical Results 


In this section, we carry out an extensive validation of all the mathe¬ 
matical models discussed in Sects. |^and[^ To this purpose, we divide 
the conducted simulations into two categories: 

• validation of the PNP model of Sect. in 3D axisymmetric 
geometries; 

• validation of the model reduction of Sect. HI 


5.1. Convergence analysis. The numerical schemes of Sect. have 
been implemented in Octave using the Octave-Forge package bim [1] for 
matrix assembly. The need to resort to radial and cylindrical coordi¬ 
nates is intrinsic in most of the geometries considered in the description 
of bio-hybrid devices, therefore we extended the package bim for these 
configurations and accurately validated the code with a broad range of 
test cases. Here we discuss the results of a convergence analysis carried 
out on the two dimensional advection-diffusion problem (14b), solved 
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on the square domain Qas = [1,2] x [0,1] in the r-z plane (the symme¬ 
try axis is r = 0). We consider the case with /i = 1, b = V'lp = [1,1]"^, 
c = 1, and where / and the boundary conditions enforced on are 
chosen in such a way that the exact solution is 

u (r, z) = z^ Inr. 



mesh size h. 




Fig. W reports the values of the L°°, and iF^-norm of the error 
M — as a function of the mesh size h. It is remarkable to notice that 
the numerical solution Uh shows the same convergence properties in the 
energy norm proved for the EAFE method in Cartesian coordinates 
in inmniEi]. Moreover, results clearly indicate superconvergence of 
the scheme in the and L°° norms, as with usual piecewise linear 
hnite elements (see [M], Chapt. 6). 


5.2. Voltage-clamp stimulation: validation of the domain re¬ 
duction. As a hrst numerical experiment, we consider a conhguration 
comprising the entire electrolyte bath surrounding the cell, in order to 
demonstrate that the main phenomena occur in the electrolyte cleft 
chosen as computational domain in Sect. The geometry, represented 
in Fig. |ll(ay is a cross section in the r-z plane of the three-dimensional 
computational domain. Boundary conditions are enforced according to 


the framework of Sect. |2.3| for the bath and for the coupling conditions, 
while on F^j^ homogeneous Neumann conditions for both potential and 
concentrations are considered, as required by the axial symmetry. The 
mesh used in the numerical computations is shown in Fig. 11(b) and 


is characterized by a local rehnement in the regions close to the cell 
membrane. In the cleft between the cell and the chip (thickness 6j = 
100 nm) we use a structured mesh in order to independently control the 
mesh characteristic dimension in the r and z directions. This allows 
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(a) Geometry 



(b) Mesh 


Figure 11 - Left: 2D geometry of the electrolyte surrounding the 
cell (cross section in the r-z plane). Dimensions: R =10pm and 
X = R. The cell is approximated as an ellipsoid with the major 
semiaxis equal to R and the other one equal to R/2. The cleft is 
the line between Tceii and Tgufe: its height is Sj = 100 nm. Right: 
computational mesh, refined all around the cell (in the zoom of a 
part of the cleft zone: the mesh is structured and refined at the 
boundaries). 

the use of very stretched triangles to achieve a more detailed descrip¬ 
tion of this area as required by the geometrical multiscale nature of the 
problem in which the ratio between cell radius and 6j is 10^. 


Parameter 

Symbol 

Value 

Intracellular potassium concentration 

„int 

140 mM 

Intracellular sodium concentration 

Ant 

^Na 

4 mM 

Intracellular chloride concentration 

Ant 

^Cl 

144 mM 

Extracellular bath potassium concentration 

Aath 

5 mM 

Extracellular bath sodium concentration 

Aath 

^Na 

140 mM 

Extracellular bath chlorine concentration 

Aath 

^Cl 

145 mM 

Potassium conductance 

9m 

250 Sm-2 

Membrane specific capacitance 

Cm 

1 pF cm“^ 

Substrate specihc capacitance 

Cs 

0.3 pF cm“^ 

Initial transmembrane potential 

Vcell - 

-85 mV 


Table 1 - Model parameter values considered in the simulations 
involving passive HEK cells. 


We consider a voltage-clamp conhguration in which the cell is stimu¬ 
lated by varying the intracellular potential from —85 mV to 50 mV, and 
we describe the transmembrane currents with the Goldman-Hodgkin- 
Katz model illustrated in Sect. |2.3[ We consider a human embryonic 


kidney cell HEK293, as in |Hll36], expressing just potassium channels, 
and the adopted parameter values are reported in Table 

Fig. 1^ shows the obtained results in terms of potential and ion 
concentrations. As a consequence of the change of the intracellular 
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Figure 12 - Spatial distribution of the electric potential and 
of the ion concentrations Ci, at the end of the transient resulting 
from a voltage-clamp depolarization at Vceii = 50 mV. 


potential, the K"*" channels open and K"*" ions flow in the extracellnlar 
space. The variation of the potassium concentration is not very large 
outside the cleft region because ions can diffuse to the bulk of the solu¬ 
tion and electroneutrality is restored in a few nanometers. In the cleft 
region, instead, ion diffusion is limited by the presence of the substrate 
and ions are forced to follow a radial pathway towards the bulk of the 
solution. As a consequence, a much higher potassium concentration is 
determined and consequently chloride and sodium ions are respectively 
attracted and repelled from the cleft region. The determined electric 
potential variation shows a parabolic profile with a peak value of ap¬ 
proximately 3 mV at the center of the cell. In view of these results, we 
can therefore state that the approximations introduced in Sect. 2.1 are 


sound. Fig. [^reports the results obtained with a simulation conducted 
only underneath the cell where it is possible to spot the steep layers at 
the membrane due to capacitive coupling and charge screening effects. 


Having started our analysis from a single cell on an electronic sub¬ 
strate, we now investigate configurations with more than one cell and/or 
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(a) tf 



(b) cci 



(c) CK 



Figure 13 - Spatial distribution of (/? and q in the portion of 
electrolyte between cell and substrate, at the end of the transient 
resulting from a voltage-clamp depolarization at Vceii = 50 mV. 


more than one electrode, focnsing on the mntnal inflnence between 
neigbouring devices. 


5.3. Voltage-clamp stimulation: signal on a neighboring elec¬ 
trode. In an ideal device for biological signal recording only the closest 
electrode to the excited cells should be recording a signal, but this ac¬ 
tually does not occur in realistic experiments. In order to investigate 
this non ideality effect, we consider a conhguration in which a cell is 
stimulated and the resulting signal is probed by two electrodes, one just 
below the cell and the other facing the bath electrolyte at a variable 
distance W. The representation of the computational domain with the 
two electrodes Tsi and rs 2 and the overlying cell is shown in Fig. 14 


The boundary conditions on T \ Te/ are the same as in Sect. 5.2 and 
the parameter values are reported in Table [T] To account for the whole 
electrolyte surrounding the cell, we enforce on Tg/ the Robin bound¬ 


ary condition described in the modeling procedure of Sect. 2.3, setting 


C* = b ■ 10“^F m“^ and v* = 10“^ms“^ in all our computations. 

The output of the simulations are the potentials measured by the two 
electrodes as a consequence of the voltage-clamp stimulation. Fig. 
reports the computed potential prohle in the cases W = 2 pm and 
W = 10 pm. Far from the cell adhesion region r > 10 pm, the po¬ 
tential decays to the bulk value ((^ = 0) with a trend proportional to 
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Figure 14 - Computational domain with two electrodes for a 
voltage-clamp stimulation recording. Figure not in scale: 6j = 
lOOnm, iicez; = 10 pm, = 1.5 pm and I = 2 pm. 



(a) W = 2 pm (b) W = 10 pm 


Figure 15 - Spatial distribution of the potential ip with different 
domains: the distance between the cell and the second gate is on 
the left W = 2 pm and on the right W = 10 pm. In both cases: 
Sj = 100 nm. 


1/r. In the first device configuration the second electrode is located 
in the region with r G [12,15] Rm and the probed electric potential 
is signihcantly different from the bulk value. In the second conhgu- 
ration, instead, the electrode is placed between [20, 23] pm, where the 
perturbation is almost equilibrated, so the measured signal is very low. 

The computed concentration prohles are shown in Fig. 16, and no¬ 


tably the prohles are in very good agreement with the results of Sect. 5.2 


Since the potential (p in proximity of the second electrode does not 
attain a uniform value, we introduce the local average 

AKs = - ''^2) dq 

d s2| Jrs2 

to be interpreted as the value probed by the second electrode and re¬ 
turned as output. We consider values of W in the range between 2 and 
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Figure 16 - Spatial distribution of the computed concentrations 
Ci with W = 2 pm. 



Figure 17 - Probed voltage AVs 2 as a function of W. Results 
for three different values of the cleft thickness Sj, obtained with a 
depolarizing pulse with Vceii = 50 mV. 


20 Rm and three different values of the cleft thickness dj, namely 50, 
100 and 150 nm, In Fig. 17 we report the computed values of AVs 2 and 
we observe an almost exponential decrease of the signal when consider¬ 
ing an increasingly distant electrode. Moreover, we observe in Fig. 


that in configurations characterized by a smaller value of 6j, a more in¬ 
tense variation of the potential is registered in the portion of electrolyte 
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Figure 18 - Spatial distribution of the potential ip with different 
domains: the distance between the cell and the second gate is set 
at VF = 2 pm, but on the left we have a cleft thickness 5j = 50 nm 
and on the right 5j = 150 nm. 


under the cell, due to the fact that higher values of the potassium con¬ 
centration occur there. From Fig. [T^we also notice that the value of 
AV's 2 decreases with the cleft thickness. A physical explanation for this 
latter result can be provided by resorting to the dehnition of electrical 
resistance for the cleft, which applies here since the electrolyte solution 
is an electrical conductor. We have 

D _ 

^el Pel~p 5 
Oel 

where pei is the electrolyte resistivity, Lei and Sei are the length and the 
cross sectional area of the cleft, this latter being linearly proportional 
to the thickness Sj. A smaller value of 6j hence results in a larger 
resistance Rei, determining a more pronounced decay of the potential 


along the cleft radius, as shown in Fig. 18(a) 


5.4. Voltage-clamp stimulation: effect on a neighboring cell. 

When more cells are attached to a bioelectronic device, if a stimulation 
is applied to just one cell, the perturbation can be sensed also by the 
neighboring ones. Here we investigate this effect by studying the device 


conhguration of Fig. 19, where a second cell on the right (denoted from 


now on as cell B) is located at a distance W from the stimulated one 
on the left (denoted from now on as cell A). 

The scheme of Fig. is not characterized by any rotational symme¬ 
try, so in principle it cannot be described using the geometrical frame¬ 
work introduced in Sect. and adopted in the previously presented nu¬ 
merical results. Our approach consists of dividing the computational 


domain of Fig. 19 into two parts along the artihcial interface Tinterf 
and assuming axial symmetry to hold separately on each of the two 
subdomains. This allows us to formulate the mathematical problem 
in each subdomain using cylindrical coordinates as in Sect. |^so that 
the solution of the problem in the whole domain is obtained through 
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Figure 19 - Computational domain with two cells and two elec¬ 
trodes. Figure not in scale: 6j = lOOnm, = 10 pm and 
Rsub = 1.5 pm). 


subdomain coupling across Tinterf using the substructuring techniques 
described in [371 [mm [9]. Since the perturbation due to cellular stim¬ 
ulation is not symmetric with respect to the axis of cell B, the above 
described approach introduces a certain level of approximation, which 
is expected to become less signihcant as the distance between the two 
cells is increased. 

We study an electrophysiological conhguration in which the inter¬ 
nal potential of both cells is controlled and while cell A is stimulated 
applying a voltage step to 50 mV, cell B is kept at the resting value 
and the membrane current is measured. Since the conhguration differs 
from that considered in Sects. |5.2| and 5.3, the resting voltage Vceq at 
which the cells are in equilibrium, (i.e., when the overall current howing 
through the membrane is zero), is not expected to coincide with the 
value of —85 mV attained in the case of the single cell conhguration. 
The new value of Weo can be determined by solving the PNP system 


subject to the condition Wi = W 2 = Weg, 
quantity, and to the following constraints 

1 f . 1 


Vceq being an unknown 


cl 




'r,2 


/”^(W2,^|r.2)^7 = 0. 


|rc2| 

The obtained value for Vceq is about —90 mV, slightly more negative 
than the previous value of —85 mV. This is probably due to the fact 
that in the case of the two-cell conhguration, potassium concentration 
in the cleft region is, on average, larger than in the case of the sole 
cell A, because of potassium current injection also from cell B. As a 
consequence, a larger dihusive hux tends to drive positive ions from the 
electrolyte cleft towards the intracellular sites so that a more negative 
clamp voltage is needed to increase the transmembrane electric held 
required to counterbalance the increased dihusive hux and restore the 
equilibrium condition of zero transmembrane current how. 

shows the spatial distributions of the potential ip and of 


Fig. 20 


the potassium concentration for two diherent values of the distance W 
between the two cells. The channels of cell A are always open and 
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(a) with W = 20 jim (b) cx with = 20 jim 



Figure 20 - Spatial distributions of ip and of cx at the end of 
the transient resulting from a voltage clamp stimulation of cell A 
to Vceii = 50 pV, keeping cell B at Vceq- Results for two different 
distances between the cells: W = 20 pm and W = 70 pm. 


by Fig. 20(a 
cell B. In Fig. |20(b) 


are injecting a current in the electrolyte, because of depolarization. 
This causes an increase of K"*" and of ip in the considered domain, which 
may lead, in turn, to the opening of the channels of cell B. In the case 
where the two cells are close enough (at a distance W = 20 pm) the 
electric poten tial exhibits a signihcant boundary layer at Fc 2 as shown 
to which corresponds an opening of the channels of 
we observe an evident depletion in the spatial 
distribution ck under Fc 2 . This is due to the fact that the potassium 
current here is entering into cell B: as physically expected the potassium 
is injected by one cell (cell A) and collected from the other one (cell B), 
in a manner that resembles, using an electronics analogy, the working 
principle of a solid-state transistor (see the series of articles HSIESIIII]). 
In the case of a larger W (Figs. |20(c) and 20(d)), the value of the 
potential in the electrolyte is lower and there is practically no current 
entering into cell B, because in this conhguration ions are free to flow 
in a larger portion of electrolyte. 

In Fig. ^ we report the integral averages of the transmembrane 
current densities and flowing out of the hrst cell and into the 
second cell, respectively, as functions of the distance W between the 
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two cells. While is unaffected by W, exponentially decays with 
W. 



Figure 21 - Computed transmembrane current densities j*'", i = 
1, 2, for the two cells as a function of W. With j*™ we denote the 
outward current flowing from the interior of cell 1 to the electrolyte, 
while with we denote the inward current flowing into the second 
cell. Results obtained with a depolarizing step at Ki = 50 mV and 
keeping 1^2 = Fceg- 


5.5. Cells with active channels. In all the numerical experiments 
described so far, we have dealt with cells with passive channels. We 
now consider the case of a cardiac cell as the one studied in [32|, which 
expresses voltage-gated ion channels. In order to describe the behavior 
of such kind of cells we adopt the well known Hodgkin-Huxley model 
(see for the details [211231129]), to reproduce the opening/closing dy¬ 
namics of the channels, so that the transmembrane ion currents are 
computed at every time t as: 

(jK = 9Kn^ {Vceii - - yK\r^J 

(15) < ]Na = QNahm^ {Vcell “ “ '^Ya|r,,J 

W/eafc 9ci (ycell ^Ircei! ’ 

where Vk, Vva and Vci are the Nernst potentials of potassium, sodium 
and chloride ions at Tceii, respectively, while and are the max¬ 
imum values of potassium and sodium ion conductances. The value 
represents the probability of a K"*" channel to be open, while the 
probability that a Na’'" channel is open is given by the product m^h. 
The gating variables n, m and h need to be consistently computed by 
solving the following system of ODEs at each of the mesh nodes on the 
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side Tceii representing the cell membrane: 


r dn _ 

d( “ 


(16) 


n] 


Pnn 


— = 0:^(1 - m) - Prnm 

^ = ah{l -h) - /3hh. 


Here we study a benchmark patch-clamp experiment in voltage- 
clamp conhguration, in which a depolarizing pulse is applied from the 
resting state (—70 mV) holding Vceii at 15 mV and considering the pa¬ 
rameter values listed in Table [U 


Parameter 

Symbol 

Value 


Intracellular potassium concentration 

^int 

140 mM 


Intracellular sodium concentration 

Ant 

^Na 

10 mM 


Intracellular chloride concentration 

Ant 

^Cl 

20 mM 


Extracellular bath potassium concentration 

Aath 

5 mM 


Extracellular bath sodium concentration 

Aath 

^Na 

145 mM 


Extracellular bath chloride concentration 

Aath 

^Cl 

150 mM 


Max potassium conductance 

9k 

00 

O 

o 

S m“^ 

Max sodium conductance 

9Na 

600 ■ 10^ 

S m“^ 

Chloride conductance 

9ci 

1.5 ■ 10^ 

S m“^ 

Membrane specific capacitance 

Cm 

1 pF cm' 

-2 

Substrate specihc capacitance 

Cs 

0.3 pF cm ^ 

Initial transmembrane potential 

Vcell - 

-70 mV 



Table 2 - Model parameter values considered in the simulations 
involving active cardiac cell. 


Fig. 22 shows the computed opening probability profiles of and 
Na^ channels, respectively, 10 ms after the onset of the stimulation, 
and as expected we can observe that the potassium channels are open 
while the sodium channels are already inactivated. The limited spatial 
variability of the obtained prohles agrees with the fact that in the 
considered patch clamp conhguration, the potential is modihed in the 
whole cell domain, so the channels on the membrane respond almost 
uniformly. Fig. ^ shows the spatial distributions of potential and 
ion concentrations of the species with active channels (K’*' and Na’'"). 
Notice the occurrence of sharp boundary layers at the cell-electrolyte 
interface due to sodium channel inactivation. 

In order to compare our results with the voltage clamp experimental 
measurements on cells with active channels, we compute the integral 
average of the channel opening probabilities and hm? over Tceii and 
we show their evolution in Fig. 24(a)[ These results are in excellent 
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(a) n^{r) (b) hm^{r) 


Figure 22 - Computed spatial distributions of the K''" and Na’*’ 
channel opening probabilities for every point of the boundary Tceii 
representing the cell membrane, 10 ms after the onset of a depo¬ 
larizing voltage-clamp pulse at Vceii = 15 mV. 





Figure 23 - Spatial distributions of the potential and of the potas¬ 
sium and sodium concentrations in the electrolyte cleft under the 
cell. Results obtained after 10 ms of a depolarizing pulse keeping 
Vceii =15 mV. 


agreement with the behavior characteristic of the classic HH model 
(see. e.g., [151 ESI EH]), but, while this latter is usually applied in a 
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lumped manner, our model has the feature of describing possible spa¬ 
tial inhomogeneities in the activation (cf. Fig. 23). In Fig. |24(b) we 
report the average ion channel transmembrane current densities com¬ 
puted according to (15), and, as expected, the typical current prohle 
after a voltage-clamp depolarization is recovered. After the onset of 
the voltage step, an inward current is induced by the opening of Na^ 
channels, which are eventually inactivated, and for longer times the 
channels open and determine a stationary outward flux. 




(b) Transmembrane currents 


Figure 24 - Time evolution of the gating variables and ion chan¬ 
nel opening probabilities (left), and of the ion current densities 
(right). Results averaged over Tceii and obtained with a depolar¬ 
izing voltage pulse at Vceii = 15 mV. 




(b) QNa = 


Figure 25 - Temporal variation of the integral mean over the 
boundary Tc,eu of the conductances of the potassium and of the 
sodium channels, computed as in (15). Results obtained with four 
different depolarizing pulses keeping Vceii = -30,-15, 0, -1-15 mV. 


As a hnal result, we consider the effects of different depolarizations 
of the cell over the ion transmembrane currents. In Fig. 2^ we show 
the variation in time of potassium and sodium conductances in corre¬ 
spondence of four different values of Vceii- We observe that gNa turns 
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on more rapidly than qk- Moreover, the Na+ channels begin to close 
before depolarization is turned off, whereas the K"*" channels remain 
open as long as the membrane is depolarized. As in the classic HH 
theory, when the cell is depolarized the Na"*" channels switch from the 
resting (closed) to the activated (open) state and then, if depolarization 
is maintained, the channel switches to the inactivated state again. 

5.6. Reduced order models. In this section we illustrate the results 
obtained solving the reduced models introduced in Sect. 

5.6.1. Validation of the 2.5D model against the 3D PNP model. In this 
section we verify the accuracy of the reduced model of Sect. |3.1| in 


the study of the same biophysical conhguration analyzed in Sect. p?2 
with the 3D axisymmetric PNP model. Since the geometrical setting is 
axisymmetric, we can reduce ourselves to a one dimensional manifold 
Dradi describing the variation of the quantities of interest along the 
sole radial direction. The considered domain Vt^ad is the union of two 
different parts U 03/ (the part where the cell is attached and the 
free part of extracellular fluid). For the approximation of the boundary 
layers introduced in Sect. m based on the simulations of the previous 
sections and on asymptotic analysis (see [32]), we set H ~ 2XDebye — 
1.6 nm. 

Simulation results obtained with the reduced model are shown in 
Fig. 26, where we immediately notice that the capacitive couplings 
obtained with 3D simulations of Sect. 5.2 are well reproduced. As ob¬ 
served in Sect. |3T and in the spatial distributions of Sect. |^, the 


quantities of interest in 0/ are not expected to sensibly vary along 
This solution behavior is conhrmed by the radial distributions of 
potential ip and concentrations q are perfectly superimposed 


z. 


Fig. 26 


on the distributions of the top and bottom quantities in the free part. 
We observe a very fast decay in 0/, as expected and as the one ob¬ 
tained in the simulation results shown in Fig. 

Having computed the averaged and the top and bottom values of the 
dependent variables, we can also reconstruct their ^-dependence, at a 
hxed point f, using the following post-processing formulas: 

(17a) 


(p{r,z) = (piz,r) + 
(17b) 

Ci{f,z) = Cjexp ( -2;* 


Ptop{r,z) - (p{r,z) 


H 


+ 




(pir,z) - (pbot{r,z) 

H 


(p{r,z) -(p{r,z) 


V 


th 


where Di and O^ are the boundary layer subdomains (see Fig. [^. The 
reconstruction of (p{r = 0,z) and of ck/ = 0,z) is illustrated in Fig. 
27 and we see that the major term in the electrolyte system behavior 
is the capacitive coupling with the cell membrane. 
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(c) K (d) Na 


Figure 26 - Spatial distribution of Tp and of Ci in the domain 
^rad = U (results obtained with = 2 • = 

2-Rceii =20pm). To account for the boundary layers: distributions 
of the top and the bottom values ptop, Phot, c°^ and c^°*. 




(a) ip{r = 0,z) 


(b) CK{r = 0,z) 


Figure 27 - Distributions along the 2 -direction of potential and 
potassium concentration at r = 0 using 


The above obtained results allow us to conclude that the reduced 
mathematical model of Sect. 3T is valid, because its predictions fa¬ 
vorably agree with those of the 3D model, but with a much smaller 
amount of degrees of freedom. 
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5.6.2. Validation of the 2.5D model against exact solutions. In this sec¬ 
tion we analyze a biophysical setting similar to that presented in |36j . 
where the authors derive the exact prohles of potential and ion concen¬ 
trations in radial coordinates under the assumption of axial symmetry. 
Their model refers to the middle plane of a cleft between a cell and a 
substrate as in the general setup described in Sect. but under the 
following simplihcations: 

• ions are assumed to flow only in the radial direction, so that 

= ji,z = 0, and no spatial dependence on cj) and 2 ; is assumed. 
The radial coordinate r is taken in the interval [0, i?] {R being 
the cell radius); 

• the influx of K’*' ion charge per volume and per time is given 
by = fx'/^ji where 5j is the cleft width and iff" is the 
potassium current density through the membrane, here assumed 
to be constant. Therefore in Eq. (6a) we only have = Xk] 

• there is no influence on the flux of ions inside the cleft from 

the two interfaces (the cell-cleft and the chip-cleft interfaces), 
neglecting the capacitive couplings described in Sect. |2.3[ This 
gives = 0 in Eq. (6c); 

• only the stationary case is considered, meaning that all quanti¬ 
ties are independent of time. 




(a) Potential (b) Concentrations 


Figure 28 - On the left: radial profile of the potential y (r). On 
the right: radial profile of the changes of ion concentrations with 
respect to their bath values Ci (r) — Results obtained with 

a source term \k = llpApm“^ and a cell radius R = 15 pm as 
in [36] . 


Fig. 28 shows the computed distributions of potential and ion con¬ 
centrations under the cell. By inspection on the analytical solutions 
reported in [SB] it can be seen that our results are in excellent agreement 
with the latter solutions. Convergence of the hnite element solution as 
a function of the radial mesh size h is reported in Fig. ^ where a 
quadratic rate can be observed before the occurrence of error satura¬ 
tion due to the nonlinear solver tolerance. This result, obtained in the 
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solution of a nonlinear problem, confirms the convergence analysis of 
Sect. 5.1, carried out on a linear model problem. 



Mesh size h 


Figure 29 - Log-log plot of the maximum norm of the discretiza¬ 
tion error ip — (ph- 

Under the hypotheses illustrated above, the potential variation at 
r = 0 is quite small (less than ImV) and the absolute changes of ion 
concentrations for CU and Na’'" are quite small too, except for K"*" ion 
concentration, from 5 mM to 8 mM. 




(a) Potential p{t,r = 0) (b) Concentrations Ci{t,r = 0) — 

„bath 


Figure 30 - Time variation of the potential and of the concen¬ 
trations with respect to their bath values, at the center of the 
junction. 


We also conduct a time dependent simulation of this experimen¬ 
tal setup, redehning the transmembrane current as Xk (t) = XrH (t), 
where Xk is the constant current used in the static simulation and 
H (t) is the Heaviside function. With this mathematical dehnition of 
potassium injection, we consider an instantaneous opening of the 
channels at t = 0, which leads to a time variation of the quantities 
and Ci (see Fig. 30, where we show the variation of these functions evalu¬ 


ated at r = 0). Transients are exhausted in about 150 ms, in agreement 
with the results of |15], while the steady-state values of ip and q agree 


well with those computed in the static case shown in Fig. 28 
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5.6.3. Validation of the Area-Contact model. In this concluding section 
we compare the results obtained with the Area-Contact (A-C) model 
proposed in Sect. 3.2 with the results of [8], which we refer to for all 
physical data and details of the electrical equivalent circuits used to 
determine the time evolution of ion concentrations in the cleft. To this 
purpose, we conduct a hrst simulation considering given concentrations 
Ci constant in time and only solving Eq. (13) (2D electrical model). 
Then, we conduct a second simulation accounting for ion dynamics, 
by adding an ODE system for the ionic concentrations (2D electrod¬ 
iffusion model). In both cases the integral mean Vj is computed as 
in (13c). The considered electrophysiological experiment is a voltage 
clamp stimulation with the depolarizing pulse shown by Fig. |31(a) and 
the values of model parameters are the same as in jH]. 



(a) Intracellular potential Vceii 
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(c) Vj electrodiffusion model 


(d) Nernst potentials VJq 


Figure 31 - (a): depolarizing pulse of intracellular potential Vceii- 
(b): integral mean of the cleft potential y(a;, y) obtained with the 
electrical model, (c): integral mean of the cleft potential (p{x,y) 
obtained with the electrodiffusion model, (d): changes of the 
Nernst potentials Vjq between junction and bath. 


As demonstrated by Fig. 31(b), the 2D electrical model accounts 


only for the fast response of the system, with a dynamics determined 
by the electrical time constant r = {Cm + Cs) /cr ~ 1.0944 pS (cr = 
qJZiLi being the global cleft conductance): when the cell 

is depolarized, almost instantaneously the potential goes to a value 
around 2 mV. The 2D electrodiffusion model, instead, describing the 
time variation of Cj, accounts also for the slow component, as shown 
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(a) ip electrical model (b) ip electrodiffusion model 


Figure 32 - Spatial distribution of the cleft potential ip in the 
circular domain at t = 0.35 ms (cell just depolarized: on the left, 
electrical model, on the right, electrodiffusion model 


by Figs. 31(c) and 31(d) both potential and concentrations have tran¬ 
sients with a time constant in the order of milliseconds, as expected 
(we have expressed the changes of extracellular ion concentrations in 
the junction as Nernst potentials between junction and bath). The in¬ 
tegral mean of the electrical potential ip increases fast to a value around 
2.5 mV and subsequently decays to a stationary level around 1.5 mV 
and the potassium concentration increases from 5 mM to 17 mM, giv¬ 
ing a Nernst potential Vf ~ —27 mV in the junction. Notably, all 
results are in excellent agreement with those reported by Brittinger 
and Fromherz in [8]. Finally, the spatial distributions of the potential 
ip computed by the 2D electrical and electrodiffusion models, and used 
to determine Vj in (13c), are reported in Fig. 32 The resulting para¬ 


bolic shape is in very good agreement with the behavior shown by the 
3D results of Sect. 13.21 


6. Conclusions and Future Perspectives 

In this article we have addressed the mathematical modeling and 
numerical simulation of ion electrodiffusion in bio-hybrid devices. This 
subject is of paramount importance in the wider scientihc context of 
neuroelectronics, where the main aim is to actually realize devices con¬ 
sisting of the integration of biological tissues with solid-state integrated 
electronic circuits. 

In this treatise we have illustrated a suitable mathematical charac¬ 
terization of bio-electronic interfaces, investigating different possible 
modeling hypotheses on the coupling between the two different envi¬ 
ronments (cell and electronic device) and on the derivation of model 
dimensional reductions, performed to decrease the computational sim¬ 
ulation effort. A hierarchy of multiscale models has been therefore 
presented and extensively validated with a broad range of numerical 
computations, obtaining sensible results and comparing them with lit¬ 
erature and experiments. 
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This mathematical description has also been applied to complex con¬ 
figurations and has proved to be able to simulate the interactions be¬ 
tween multiple cells and multiple devices. Even if the present work is 
not a faithful copy of a real-world biophysical setting, it can be con¬ 
sidered a first step for the construction of mathematical models to be 
used in the design of actual devices. 

Clearly, future research is needed to provide a better description of 
the complex multiscale/multiphysics problem object of our investiga¬ 
tion. Among possible developments, we mention: 

• a more accurate modeling of the electronic substrate, which can 
be useful in studying different types of stimulation, for example 
a different polarization of the chip influencing the cell; 

• a model for the chemical binding mechanism of the ions to the 
electronic substrate is also required, in order to fully describe 
the EOSFET device; 

• a coupling between electro-chemical and fluid-mechanical sys¬ 
tems, in order to account for the forces due to pressure differ¬ 
ences and flow in the aqueous medium; 

• a more realistic description of the problem geometry, with full 
three-dimensional computations including the intracellular fluid 
can be useful to faithfully reproduce the entire phenomena; 

• an application of the computational model to the simulation of 
electrophysiological experiments in current-clamp conditions. 

The above mentioned improvements, particularly, the study of the 
current-clamp protocol, should give the realistic chance to go further 
in the study of the interactions between multiple cells, maybe intro¬ 
ducing a neural network and simulating a whole brain slice, as in the 
experimental results of [26l 07] . 
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